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We investigate the dynamical correlation function of a quantum-mechanical two-state system which 
is coupled to a bosonic heat bath, utilizing the equivalence between the spin-boson Hamiltonian 
and the 1/r 2 Ising model. The imaginary-time correlation function is obtained from Monte Carlo 
simulations on the Ising system and then continued to real time by a Pade approximation. In the 
unbiased system, the transition from oscillatory to strongly damped behavior is found to occur at 
a coupling strength close to a = 1/2. The biased system favors coherent relaxation and displays a 
significantly larger crossover value a c . We introduce the quasiparticle picture to describe the relevant 
behavior at intermediate time scales. Within this approximation, we map out phase diagrams for 
the unbiased and biased systems. 
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I. INTRODUCTION 

The quantum mechanical two-state system, coupled to 
a dissipative environment, provides a universal model 
for many physical systems. An important example of 
current interest is the problem of defect-tunneling in 
solids .cl Our interest in this problem was renewed when 
connections between current topics in high temperature 
superconductivitycl and quantum Hall efTectQ were per- 
ceived. In fact, all systems that can be described by a 
double-well potential associated with a generalized coor- 
dinate, with appropriate restrictions on the parameters, 
reduce to a two-state system at sufficiently low tempera- 
tures, when only the ground states of the two individual 
wells remain relevant. □ 

At low temperatures, we can choose the localized 
ground states of the two potential wells, in the absence 
of tunneling, as the basis of our two-dimensional Hilbert 
space. The overlap of the two wavefunctions leads to 
quantum mechanical tunneling, described by a transition 
matrix element between the two effective states of the 
system. Decoupled from its environment, we can write 
the Hamiltonian of the system in pseudospin language as 



-Wo = ~~^ a x + 2 " 2 



(1) 



a x and a z are the usual Pauli matrices, so that A/ 2 is the 
tunneling matrix element, and e describes the bias of the 
system, i.e., the difference in the ground state energies of 
the two localized states. The isolated system is trivial to 
diagonalize and exhibits coherent tunneling between the 
two states. n 
If the system is coupled to a dissipative environmentE, 
we can expect the tunneling to lose its phase coherence. 
This can happen even at zero temperature if the con- 
tinuous spectrum of the macroscopic dissipative environ- 
ment extends down to zero frequency. A particularly 
elegant model of the environment has been known for 



some timeEl In this model, a set of harmonic oscillators 
are coupled linearly to a z . The full Hamiltonian is 

A e + 



H 



2 2 
^-^2fa(ai +a a ), 



(2) 



where the a Q are destruction operators of the harmonic 
oscillators with frequencies u a . The quantity f a repre- 
sents the coupling strength of the two-state system to the 
coordinate of the a oscillator. This is also known as the 
spin-boson HamiltonianQ that was proposed in the con- 
text of a model of a local magnetic impurity coupled by 
spin- flip scattering to the conduction electrons of the host 
metal, known as the Kondo problem. The low-energy 
particle-hole excitations of the conduction electrons de- 
fine the oscillator states. The environment is completely 
characterized by its spectral density J(to) 



(3) 



As a model of a linear dissipative environment, we con- 
sider an Ohmic environment that classically exhibits fric- 
tion of the fbrm F — —rjq, and is described by the spec- 
tral densitya: 



J ^ = { 0, c»c c 



(4) 



The assumption of an Ohmic bath reduces the charac- 
teristics of the environment to two parameters, the di- 
mensionless coupling strength a and the characteristic 
cutoff scale uj c . The generality of this model is not ob- 
vious, unless the phenomena to be studied are strongly 
dominated by low energy processes. Indeed, the results 
of the present paper can be trusted for only such low en- 
ergy processes, which are akin to dynamic quantum crit- 
ical phenomena. We shall be interested in the dynamical 
correlation function 
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C(t) = -<[<*(<), <r,(0)]>, 



(5) 



or its Fourier transform, the response function ^"(o;), at 
zero temperature. In particular, we will consider its qual- 
itative and quantitative dependencies on the parameters 
A, e and a, appearing in the guise of dimensionless ratios 
A/lo c and e/to c (h and kg are set to unity). 

In an analogous model, one can also couple a two-state 
system to a Fermionic environment that can be described 
by a non- interacting set of quasiparticles,u or to a Fermi 
liquid. Except for certain restrictions on parameters, this 
is equivalent to the oscillator model for the low-energy 
states of the tunneling degree of freedom. More intrigu- 
ing, however, is the case where the tunneling degree of 
freedom is coupled to a Luttingcr liquid as in a quasi- 
one-dimensional election gas, or to the edge states in a 
quantum Hall systemja but this topic is outside the scope 
of the present paper. _ 

It is known for some timeQ that the partition function 
of this system can be transformed into the partition func- 
tion of an Ising model with long-range interactions. A 
particular Ising site corresponds to the pseudospin state 
at a particular point on the imaginary-time axis. This 
connection enables us to determine the imaginary-time 
correlation function of the two-state system from Monte 
Carlo simulations of the Ising model. A Wick rotation 
then leads us from imaginary time to the above men- 
tioned real-time correlation function C(t). This proce- 
dure will be described in detail in the following sections 
and constitutes a reliable nonperturbative evaluation of 
the dynamics of the system in regimes in which dilute 
blip approximation employed previously is not reliable. 
The present paper is based on ideas introduced,in a short 
communication by Chakravarty and Rudnick.tj 

We know that the correlation function C(t) will exhibit 
coherent oscillations for zero coupling. Weak coupling 
to the environment leads to damped oscillatory behav- 
ior. For small a and A/oj c , scaling argumentsa yield the 
renormalized tunneling frequency 



A,. = A (A/u 



\a/l- 



(6) 



If the influence of the environment is strong enough, we 
expect a completely incoherent decay of the time corre- 
lation function; further increase of the coupling leads to 
broken sjaametry and the tunneling degree of freedom is 
localized.!!!] One aim of the present paper is to map out 
the coherent and incoherent regimes in a phase diagram, 
at zero temperature, and to provide quantitative results 
for the behavior near the crossover value of a. We will 
also obtain a qualitative picture of the behavior of C(t) 
in a biased system (e 7^ 0). 

In Sec. we will derive the precise correspondence be- 
tween the spin-boson Hamiltonian and a classical Ising 
spin system. Section III deals with aspects of the Monte 



Carlo simulations on the Ising model, and with the Pade 
approximation method used to obtain real-time results. 



Those results are presented in Sec. [V. Section [v] in- 
troduces the quasiparticle picture as an elegant way to 
model the essential physics at intermediate times. Re- 
sults obtained within the quas iparticle approximation are 
presented in Sees. [Vl] and |VII| . The last section compares 
our results with some exact results for the zero-frequency 
limit of the spectral function and for the Toulouse case 
a = 1/2. 



II. THE DISSIPATIVE TWO-STATE SYSTEM, 
THE COULOMB GAS MODEL AND THE 
INVERSE-SQUARE ISING MODEL 

The partition function corresponding to rthe Hamilto- 
nian (|^) can be cast into the Coulomb-gasciJ form, that 
is, a one-dimensional system of alternating positive and 
negative charges. Expanding the partition function in 
terms of a x in imaginary time, we get: 



Z = Tr 



Tr 



-0H C 



71=0 



x / dri ■ 
'0 



dr n Ha(ti) • ■■H A {r ri 



where 



A 

Ha = — ^ffi, 



(7) 



(8) 



Hoo = ^2 W a a l a a + y ( M a a + a ") + 6 ] > ( 9 ) 

a \ a / 

and Ha(t) = e rHa ° HAe~ TH °° is the operator in the in- 
teraction representation. After integrating out the envi- 
ronmental degrees of freedom, we arrive at the following 
expression for Z: 



00 /A\ 2n fP 



n=0 



dT 2n 



(10) 



where 



x j^ e -Wc "l s ~ s 'l + 2n a coshuj a (s — s') 



(11) 



m = { 



+ 1, T 2 k < S < T 2 k+l 
-1, T 2 fe_l < S <T 2 k 



2£(-l)"0( S -T„), 

71 = 



(12) 
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The algebraic manipulations are detailed in Ref. ^|. n a = 
l/(e /3 "° — 1) is the Bose occupation factor, and Z is the 
partition function of the environment. At this point we 
will introduce the explicit form of the spectral density: 



J{lo) = 7r — ^a) = 2iraoj8(u)/ui c ). 



(13) 



Here, 9(x) symbolizes a generic cutoff function: 9(x) = 1 
for x — 0, 9(x) — for x ^> 1, and smooth in between. 
Then 



o 



U = -~ J dsfa) J ds'tis') 

coshTwd- \s-s'\ 
dww — 1 



(i 



sinh 4p 



(wM). (14) 



An approximate treatment of the cutoff consists in re- 
moving the cutoff from the integral in this expression, 
and reintroducing it in the form of the regularization 
condition |r — t'| > ui^ 1 in the partition function (10). 
This corresponds to a hard-sphere repulsion between the 
Coulomb gas charges. Then the integral over uj yields: 



f 

Jo 



dw ' 

7T 



cosh[w(f - \s-s'\ 
' sinh%i 



sui I —\s - s 



(15) 



Next we can perform a two-fold partial integration over 
s and s' , noting that the boundary terms are zero and 
d£/ds = 2^ fe (— l) k S(s — rfc), arriving at: 



Z = Z 



n=0 



At c 



In r j3 



T c JO 



dro 



dr-2 n 



+2a 



i<3 



expl e] 




j 


In 


A 




7TT C 



2)1 



/3 



75) 



(16) 



The r c = 1 /w c appearing in the upper boundaries of inte- 
gration is the regularization condition that replaced the 
high-frequency cutoff in J(u>). In this picture, the posi- 
tions Tj of the charged particles correspond to spin-flips 
of the two-state system on the imaginary-time axis. It 
should be emphasized that this treatment of the cutoff is 
somewhat artificial, and we therefore do not have a com- 
plete description of the high-energy details of the system. 
Note, however, that the logarithm in the above expres- 
sion (16) vanishes as r — r' — * r c , so that our treatment 
is self-consistent. A more rigorous treatment would lead 
to additional terms of the form (Tj — Ti)~ x ,x > 1 in the 
Hamiltonian. The effects of these terms on the results for 
intermediate time scales can be summarized by replacing 



At c with an effective value (Ar c ) e ff. Otherwise, the cut- 
off procedure, does not affect our results for intermediate 
time scales.l!3 Furthermore, since the In is the only scale 
invariant term in the series, the additionaL-terms cannot 
change the critical behavior of the system.lla 

Let us now consider a one-dimensional Ising model 
with long-range interactions and periodic boundary con- 
ditions, given by: 

Z= J2 <'*!>{ •' : '- S '< S '.< h zl S }- 

Si—Sn j<i j 

Following Cardyjll we can rewrite the partition function 
in terms of interactions between spin flips or kinks: 



zZv 2n 

n=0 



dr\ 
a 



dri 



dr-2 



exp{£(-l) 



-aB-tfC?)}. 



(18) 



where U (k) is defined by 

V(x) = U(k + l)-2U(k) + U(k-l), (19) 

and y = e 2U ^ is the chemical potential, or fugacity, of 
the system. 

The short-distance cutoff at t = a was introduced to 
recapture the high-energy properties of the discrete lat- 
tice, after allowing the kinks to move in the continuous 
interval r = ■ • ■ Na. Here, a is the lattice constant of 
the Ising model. Again, this approximate treatment of 
the cutoff changes the high-energy behavior of the sys- 
tem, but leaves the low-energy physics invariant, aside 
from a modification of parameters related to the cutoff. 

Comparison with ( |l6| ) yields a = t c and 



U(n) = -\n 



N . 7m 
— sin — 

7T N 



so that 



V(n) 



2 sm\Tm/N)' 



n > 1, 



n > 2, 



(20) 



(21) 



where we neglected terms of order (ir/N) 4 . V(l) is de- 
termined by U(0), which in turn depends on the fugacity 
y _ e 2U(o) _ Ar c /2. In the limit N -> oo we get: 



U(P)=V(l)--i, 



(22) 



where 7 = 0.577 ... is Eulcr's constant. In the familiar 
language of the Ising model, the Hamiltonian reads: 

fiiHi = — 7r- SiSi+i 



2 ^ 



Jlr \ -* 

9 =5 



{w/N) 2 SiSj 



2 f£sm A [*(j-i)/N] 



hY / S l , (23) 
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where Jlr — a, and J u 



V 



Jlr = -2V(1), so that 
LR . To summarize, the correspondence 



NN 



of the parameters is as follows: 



2uj c . 



= exp{-J NN - (1 + 7) J L r}; 



a = Jlr; 

— =h; 
/3u c = N. 



(24) 



Finally, the anisotropic spin-half Rondo, model can be 
cast into the same Coulomb gas form.0 In particular, 
the case a = 1/2 corresponds to the exactly solvable 
Toulouse limit of the Kondo problem, which constitutes 
an important check of our numerical results. 



III. NUMERICAL TECHNIQUES: MONTE 
CARLO SIMULATION AND PADE 
APPROXIMATION 

Due to the equivalence of both models, a calculation 
of the spin-spin correlation function Sjfli — j\) of the 
Ising model ( p3| ) provides us with the imaginary-time 
correlation function C(t) — (cr z (r)cr z (0)) of the two-state 
system. The former can be computed using standard 
Monte Carlo techniques. Since the regions of parame- 
ter space of primary interest to us lie in the vicinity of 
the Ising model phase transition, the accuracy of stan- 
dard Monte Carlo algorithms is greatly suppressed by the 
phenomenon of critical slowing down: The autocorrela- 
tion function, which measures the number of Monte Carlo 
steps necessary to obtain two statistically independent 
configurations, diverges at the phase boundary. To get 
a correct picture of the thermodynamics, an extremely 
high number of lattice sweeps becomes necessary. 

To circumvent|-this problem, we use the Swendsen- 
Wang algorithmJU where the Ising model is mapped 
onto a percolation model: In every Monte Carlo step, 
bonds are drawn with probability = 1 — e~ Jij be- 
tween each two parallel spins, and no bonds are drawn 
(Py = 0) if the spins are antiparallel. The spins do 
not necessarily have to be nearest neighbors, and 
is the sum of their nearest-neighbor and long-range in- 
teractions. Then, spins which are directly or indirectly 
connected by bonds are grouped into clusters, and each 
spin in a cluster is assigned the same spin state with 
a probability = e ±nh /(e nh + e~ nh ) for spin-up and 
spin-down, respectively. Here n is the number of spins 
contained in this particular cluster. It is straightforward 
to check that the principle of detailed balance is satis- 
fied for this algorithm, so that the distribution of spin 
states converges to the Boltzmann distribution. Our im- 
plementation of the algorithm performs the formation of 
clusters simultaneously with the assignment of bonds, so 
that no additional computation time has to be invested 
here. 



Next, the Fourier transform C(ui n ) at the Matsubara 
frequencies is obtained,. The Pade approximant method 
of Vidberg and Sereneta is used to continue this spectral 
function from the positive Matsubara frequencies onto 
the real axis. Essentially, the spectral function is ap- 
proximated by a rational function, with numerator and 
denominator polynomials of order N/2, where N is the 
number of Matsubara points used in the approximation. 
This rational function can then be trivially continued into 
the complex plane, given that the first quadrant of the 
plane is free of singularities. The imaginary part of the 
resulting function is the response function x"(w), as we 
shall see below. 

The Pade approximation, in the context of analytic 
continuation, is an ill-defined procedure in the sense that 
small errors in the input data can lead to severe dis- 
crepancies in the output data. Accuracy of the simula- 
tion data is therefore a critical issue. We used 8 x 10 6 
Monte Carlo sweeps on a lattice with TV = 256 Ising spins 
for most points in parameter space, but had to increase 
this number in the vicinity of the crossover value of a to 
2 x 10 7 sweeps. The results are essentially independent of 
N = [3lo c , so they converge well in the zero-temperature 
limit. However, simulations on a finite spin system will 
not allow us to capture the localization effectgl3 that oc- 
cur for a — > 1. 

To get a measure for the simulation errors, we calcu- 
lated the Pade approximant several times for each data 
point, using different numbers of lattice sweeps and Pade 
points. The results were filtered by checking the basic an- 
alytic structure of x"( w ) an d the behavior for |cj| — > 00. 



IV. RESULTS 

In this section we present results for the anti-sym- 
metrized correlation function C{t) = A ([<J z (t), a z (0)]} at 
zero temperature, as obtained from Monte Carlo simula- 
tions on the Ising system (^3|), and continued to real time 
by a Pade approximation. The analytic continuation is 
carried out in the complex frequency plane. Therefore 
we introduce Green's functions corresponding to the cor- 
relation functions in imaginary and real time as: 



Q{r) 



1 

•- Tr 
Z 



G(t)=--Tr 



e -tiH e \T\H -\r\H 



~/3H iH\t\ -iH\t\ 

e r e 'a z e 1 cr. 



-c(M), 



= -i(a z (\t\)a z (0)) 



(25) 



These Green's functions are bosonic in character (Q (r + 
0) = G{t),t < 0), so that their Fourier transforms are 
connected in the complex plane astJJ 



Q{— iui n ) = Re G(uj) + i Im G{u>) tanh 



(3uj 
~2~' 



(26) 



which leads to the following relation between C(r) and 
the response function x"(u/): 
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FIG. 1. The correlation function C(t) = ^ ([a z (t),a z (0)]} 
for A/u>c = 0.1 and different values of a (top). Oscillatory 
behavior can be observed up to a coupling strength a — 0.4 
(bottom). 
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FIG. 2. Oscillation frequency ujq and damping coefficient A 
corresponding to the curves shown in Fig. |l| 



die 1 ' 



tanh I 



t (*W(Q)]> 



dt e luJt Re(a z (t)a z (0)} 



ImC(— iu>). 



(27) 



Figure [T] shows some examplesEj for C(t) in the un- 
biased case (e = 0), for A/uj c = 0.1 and various values 
of a. We can clearly observe coherent, weakly damped 
oscillations of the form C(t) = Z sin(tjnt) e~ A '*' for small 
values of a. A more detailed plot (Fig. [I] bottom) reveals 
that oscillations are indeed visible up to a ~ 0.4. Past 
that value, the accuracy of our numerical data does not 
allow to resolve oscillatory behavior anymore. Figure § 
shows oscillation frequency lvq and damping coefficient A 
as obtained from the data in Fig. [ij 

In the case of a biased system (e ^ 0, see Fig. I), the os- 
cillation frequency is greatly enhanced, while the damp- 
ing coefficient remains largely unaffected by the bias. The 
oscillations are weaker, since the system dwells in the en- 
ergetically favored state most of the time. 
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of e/tJe- The biased system favors coherent oscillations. 
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FIG. 4. top: The imaginary-time correlation function 
C{ui n ) at the Matsubara frequencies, as obtained from Monte 
Carlo simulation (circles), for A/uj c = 0.2, a = 0.2, and its 
approximation by a Lorentzian (solid line). Bottom: Spectral 
function S(u>), as calculated from the Lorentz fit (solid line), 
and by Pade approximation (dashed line). 



V. ANALYTIC PROPERTIES OF THE 
SPECTRAL FUNCTION 

An example for the Fourier transform of the imaginary- 
time correlation function C(uj n ) — J Q dr e lu)nT C(t) is 
shown in Fig. |3| Note that, since C(r) is real and sym- 
metric in r, only the positive- frequency part has physical 
meaning. Analytic continuation to negative frequencies 
reveals, however, that this function is very well described 
by a Lorentzian: 



C(u> n ) 



(u> n + A) + i 



const. 



(28) 



This approximate form seems to hold extremely well in 
large parts of parameter space. Typical deviations are 
of the order of 1%, but become larger close to a = 1/2. 
However, we cannot expect this simple approximation to 
give reliable results near the crossover value of a, and it 
will not be used in the following to obtain any quanti- 
tative results. In particular, it fails to reproduce the ex- 
plicit form ([58j) of the spectral function for the Toulouse 
case a = 1/2. 

We consider the spectral function S(uj) — x"(u;)/a;. 
According to (p7[), the approximate form of S(uj) displays 
a simple four-pole structure in the complex plane: 



5(w) 



aQ- 1 



(29) 



Here, Q = u>o/\ is the Q- value of the resonance, luq is 
the oscillation frequency, and A the damping coefficient. 
Again, this simplified form of the spectral function can- 
not be expected to give reliable quantitative results and 
will be used solely for illustrative purposes. All results 
are obtained from Pade approximants to C(iv n ) based on 
up to 256 Matsubara points. 



The Quasiparticle Picture 

Applying the Pade approximation to obtain the precise 
form of S(uj) in the complex plane, we discover that its 
analytic structure is very similar to that of the Lorentzian 
approximation (p9f). It is dominated by four poles at 
ijj = ±uiq ± i\, as can be seen in Fig. |[ The physi- 
cal meaning of these poles requires some consideration: 
We know from the spectral representation of the correla- 
tion function that it has a branch cut along the real axis, 
and is otherwise free of singularities on the first Riemann 
sheet. If the spectral density associated with the branch 
cut is sharply peaked around a finite frequency uj$, it 
is legitimate to approximate the dominant part of the 
spectral density by a simple pole at ujq + i\, the "quasi- 
particle pole" , which gives rise to the oscillatory part of 
the correlation function: 

C(t) = Zsin(uj t) e~ A|t| + (incoherent part) (30) 

Alternatively, we can view the quasiparticle pole as re- 
sulting from an analytic continuation of the spectral func- 
tion past the branch cut onto the second Riemann sheet. 
This is done, for example, in Fermi-liquid theory, and 
accurately reproduces the intermediate-time behavior of 
the correlation function. The intermediate-time range- 
over which Fermi-liquid theory is valid, is defined byE£l 
1/e/c <C t < 1/7/c, where and 7fc are quasiparticle en- 
ergy and damping, respectively. In particular, neither 
the short-time decay proportional to t nor the exponen- 
tial long-time decay implied by Fermi-liquid theory re- 
flect the true behavior of the correlation function. 

In fact, the Pade approximation naturally yields this 
analytic continuation past the branch cut, since the ex- 
istence of a branch cut on the real axis cannot be in- 
ferred from the shape of C(lu) on the imaginary axis. 
The Pade approximant therefore produces precisely this 
"quasiparticle picture" of the spectral function. The fact 
that it does not reveal any other characteristics indicates 
that almost all of the spectral weight contributes to the 
quasiparticle peaks. As in Fermi- liquid theory, the quasi- 
particle picture will reproduce neither the high-energy 
behavior nor the long-time tail of the correlation func- 
tion. These are dominated by the incoherent part of the 
spectral function, which corresponds to other, more com- 
plicated singularities in the complex plane. 
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FIG. 5. The absolute value of the spectral function 
in the complex uj plane, as obtained from a Pade approximant 
using 256 Matsubara points, for A/uj c — 0.2, and a = 0.2 
(top), a = 0.35 (bottom). The analytic structure is domi- 
nated by four poles at uj = ±ojo ± *A. The insets show S{oj) 
on the positive real axis: the quasiparticle peak is visible only 
for a < 0.33. 
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FIG. 6. The journey of the quasiparticle poles in the com- 
plex w-plane 



We might also consider the symmetrized correlation 
function 

C s (t) = Re{a z (t)a z (0)) = \ (<J z (t)a z {0) + a z (0)a z (t)) , 

(31) 

which acquires an additional branch cut contribution (at 
T=0) from the tanh factor in Eqn. (|27 



cr(t) = 



dui 
2^ 



(32) 



Here, C(— to) is the analytic continuation of C(u>) to 
negative frequencies. In particular, this term leads to 
an asymptotic 1/t 2 decay. When Q < 1, this incoherent 
background contribution will become dominant at inter- 
mediate and short time scales as well, so that in this 
regime it might be difficult to observe coherent oscilla- 
tory behavior experimentally in the symmetrized correla- 
tion function C s (t). As in Fermi-liquid theory, the time 
interval over which the quasiparticle picture accurately 
describes C s (t), is bounded by 1/u>q <C t < 1/A, which 
implies Q 1. 

Note that these considerations do not apply to the anti- 
symmetrized correlation function C(t), since the contri- 
bution ( |32| ) is purely real. For this function, we expect 
deviations from the quasiparticle behavior only at short- 
time scales comparable to the inverse cutoff l/u) c - Fur- 
thermore, we have seen in Sec. [V that damped oscilla- 
tory behavior is observable at time scales up to a few 
times the inverse damping coefficient. At longer time 
scales, a crossover to algebraic decay may also be ob- 
served in the anti-symmetrized correlation function C (t) . 



VI. OSCILLATORY AND DAMPED BEHAVIOR 



The results presented here are based on the quasipar- 
ticle picture introduced in the previous section. We have 
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seen in Sec. |y| that the quasiparticle picture offers an 
accurate description of the intermediate-time dynamics 
at feast for coupling strengths up to a ~ 0.4. Close 
to a = 1/2, the spectral function may acquire a signifi- 
cant incoherent part due to more complicated singulari- 
ties in the complex plane, and the quasiparticle picture 
might not offer an accurate description of the dynamics 
anymore. Therefore, the results very close to a = 1/2 
should best be viewed as characteristics of the quasipar- 
ticle model. 

For a = (no coupling to the environment), each 
two quasiparticle poles coincide on the real axis, yield- 
ing the familiar oscillatory behavior of every quantum- 
mechanical two-state system. As we turn on the coupling 
(a > 0), the poles move away from the real axis into the 
complex regime, as sketched in Fig. |^, corresponding to 
weakly damped oscillations. For even larger values of a, 
the poles meet with the imaginary axis, giving rise to 
overdamped relaxation. 

The location of the poles in the complex plane was 
obtained numerically from the results of the Pade ap- 
proximation. Fig. ^ shows the real (ujq) and imaginary 
(A) part of their location as a function of a for different 
values of AJut c . Note that the results agree very well 
with Sec. IV. The real part ujq fits nicely to a power-law 
curve: too(a) ~ (a c — ot) v , a < a c . Since we cannot ad 
hoc assume that the power-law form holds for arbitrary 
a < a c , only points with a > 0.3 were included in the 
fits. With the exception of A/ui c = 0.5, the graphs follow 
the power-law form surprisingly well for all a < a c . The 
results do not change significantly if we limit ourselves to 
an even smaller number of data points in the vicinity of 
a c . 

Th-e Q value, Q = ujq /A, is shown in Fig. ^. The scaling 
forn 



S(u) 



(33) 



which holds for A <C oj c , implies that the Q value is a 
universal function of a, which is independent of A/w c . 
This is indeed the case for A/uj c < 0.2 and a not too 
close to the critical value. 

For values of a that are not too large, the oscillatory 
behavior of the system is visible in the spectral function 
in the form of inelastic peaks on the real axis, centered 
around w — zhwo. If we increase the damping, the peaks 
broaden and move closer to the origin. At Q = 1, when 
the half-width of the peaks starts to exceed their separa- 
tion, they can no longer be distinguished, and the spec- 
tral function now appears to be centered at u) — 0, as can 
easily be verified from the approximate form (|2^). The 
Q factor reaches unity when a ~ 1/3, in agreement with 
results obtained by a numerical renormalization jgroup 
calculationlij and an analytic form factor approachti Os- 
cillatory behavior persists beyond that point, however, 
due to the presence of the quasiparticle poles in the com- 
plex plane at nonzero values of luq . This was confirmed 
in Sec. IV without resorting to the quasiparticle picture. 



o co,/co c 
x A/co 




FIG. 7. Tunneling frequency loo (circles) and damping co- 
efficient A (X's) as functions of the coupling strength a, for 
A/uj c = 0.05,0.1,0.2 and 0.5, respectively, too is strictly zero 
where no error bars are shown. The solid lines are power-law 
fits. 



A/co c = 0.05 
■ A/co c = 0.1 
♦ A/(0 = 0.2 

A/co =0.5 




a 



FIG. 8. The Q value Q = Wo/A as a function of a, for 
A/lu c = 0.05,0.1,0.2 and 0.5, respectively. For A/uj c < 0.2 
and a not too close to a c , the data points fall onto a universal 
scaling curve. The dashed line is expression (B4). 
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FIG. 9. The phase diagram of the two-state system in the 
(A, a) plane. The four data points for a c were obtained from 
power-law fits to ojo(a); the dashed line is a guide to the 
eye only. The va lues f or A/u> c include the cutoff corrections 



estimated in sec. VIII 



The physical meaning of a = 1/3 is simply that the 
damping coefficient starts to exceed the oscillation fre- 
quency. 

Within the quasiparticle point of view the transition 
to strongly damped behavior takes place when the poles 
coincide with the imaginary axis (ujq = 0). For small 
A/oj c , the equivalence between the two-state system with 
a = 1/2 and the Toulouse limit of the Kondo model pro- 
vides good reasons to assume that this transition occurs 
at a c — 1/2. Calculations of the correlation function 
P(t), which is the conditional average (c z (t)} with t 
constraint <r z (t) = lir t < 0, support this conclusion 
An exact expressions for the Q value associated with 
P(t) is known in the limit A/uj c — ► 0: 



Q 



7T 

cot I — - 



2 1 



a 



(34) 



This expression matches the data in Fig. || almost exactly, 
for a not too close to 1/2. This is an indication that P(t) 
and C (t) possess an equivalent structure at intermediate 
times. The discrepancy near a = 1/2 may either indicate 
a breakdown of the quasiparticle picture, or may arise 
from the fact that we consider finite values of A/uj c . 

Our results indicate a crossover value a c that is slightly 
smaller than 1/2. We have to be aware that it is in 
general a difficult task to extract critical coefficients and 
critical parameter values from Monte Carlo simulation 
data. Simulation errors increase significantly near the 
transition point, and in our case the errors get magnified 
by the procedure of analytical continuation to the real 
axis. Nevertheless, according to our results, loq is strictly 



zero at least for a > 0.48. However, for A/oj c < 0.1 the 
Ising correlation function falls off so fast that the Pade 
approximation becomes increasingly unreliable. The cal- 
culations for A/u> c — 0.05 and a near its critical value 
were performed on a lattice of N = 512 Ising spins, but 
simulations with an even larger number would be neces- 
sary to confirm the accuracy of the results. These are 
very difficult to conduct due to the N 2 beha vior of the 
computation time. As we shall see in sec. VIII, cutoff cor- 
rections become significant for small A/u> c , and prevent 
us from approaching the limit A/ui c 0. With cutoff 
corrections and error bars taken into account, the results 
presented here allow tha conclusion that a c = 1/2 is the 
correct crossover valueEH in the limit A/uj c — > 0. 

The results of the power-law fits are summarized in the 
following table (y is the "critical exponent" of wn). The 
phase diagram of the system is sketched in Fig. |9j. 



A/lu c 


a c 


V 


0.05 


0.464 ±0.015 


1.11 ±0.11 


0.1 


0.460 ±0.012 


0.91 ±0.12 


0.2 


0.431 ±0.010 


0.63 ±0.09 


0.5 


0.38 ±0.04 


0.47 ±0.31 



VII. THE BIASED TWO-STATE SYSTEM 

We will now investigate the influence of a bias (e > 0) 
on the dynamical behavior of the system, i.e. one of 
the states is energetically favored. The spin-spin corre- 
lation function, as defined above, then acquires a delta 
function peak at zero frequency, due to the non- vanishing 
expectation value (a z (r)) . If we instead consider the irre- 
ducible correlation function (<t z (t)(j z (0)} — (o^) 2 , we dis- 
cover that its Fourier transform is again surprisingly well 
described by a four-pole structure. The bias shifts the lo- 
cation of the quasiparticle poles in the complex plane, but 
does not otherwise change the analytic structure of the 
dominant part. Fig. |l^ shows the oscillation frequency 
wo and damping coefficient A as a function of e/uj Cl for 
A/uj c = 0.1 and three different values of a. 

In the first graph, a = 0.2, the damping coefficient 
A remains largely unaffected by the bias. As we might 
expect from the undamped case (a = 0), where ujq = 
V A 2 + e 2 , the tunneling frequency increases linearly with 
the bias, in the limit of large e. This leads to an overall 
linear increase in the Q value. These features were al- 
ready observed in Sec. [Tv|. The graphs for a = 0.5 and 
a = 0.75 show the same qualitative behavior for large 
e/uo c . We see that a bias can induce a transition from 
strongly damped (ujq = 0) to oscillatory (ujq > 0) behav- 
ior. Since the system dwells in the lower-energy state for 
most of the time, the interaction with the environment is 
reduced, and damping is less effective in the presence of a 
bias. At zero temperature and for moderate coupling to 
the environment (0.5 < a < 1.0), the system will always 
display oscillatory behavior if the bias exceeds a critical 
value e c (e c /uj c ~ 0.01 for A/u c — 0.1, a — 0.5, and 
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FIG. 10. Tunneling frequency and damping coefficient as a 
function of the bias e for A/ui c = 0.1 and different values of 
a. The biased system favors coherent relaxation. 
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FIG. 11. The phase diagram in the (e, a) plane, for 
A/uj c = 0.1 



e c /uj c ~ 0.03 for A/w c = 0.1, a = 0.75). The resulting 
phase diagram is sketched in Fig. [ll| 

At long time scales, the dynamics are governed by the 
zero-frequency behavior of the spectral function. Fig. 
|l2| shows lim^^o S'(a|4 q as a function of the bias. The 
theoretical predictionca for large e/u> c is 



lim S (u>) 



4q-6 



(35) 



This is in good agreement with the results shown. 



VIII. SHIBA'S RELATION AND THE TOULOUSE 
LIMIT 

Shiba's relation was originally proven for the Ander- 
son model, hiit later generalized to the two-state sys- 
tem withouto and withE3 a bias. It connects the zero- 
frequency behavior of the spectral function to the static 
susceptibility as: 



lim i%) = -a Xo , 
where the static susceptibility 
uJcXo = lim iWcx'M 



(36) 



(37) 



can be directly extracted from the Monte Carlo simu- 
lations on the Ising system. Shiba's relation therefore 
constitutes another important check of our numerical ap- 
proach. 



>! X 



i<r 

10° 

10" 2 
10 4 

10 2 

10° 
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10 4 

10 2 

10° 
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lim co->o m c S(co /co c ) 



a = 0.2 




a = 0.5 



a = 0.75 



10" 



10" 



e/co 



FIG. 12. Zero-frequency limit of S(td) as a function of the 
bias. The solid lines (with error bars) are our simulation 
results. The 'X' show S(ui — > 0) as given by Shiba's relation. 
The dashed line in the graph for a = 0.5 is the exact analytical 
expression (^) for the Toulouse case. 
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a A/co° = 0.1 
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o V(») = 0.5 




co/A T 



FIG. 13. The spectral function in the Toulouse limit 
a = 1/2. The solid line shows the analytic expression, while 
the symbols represent numerical results for different values of 

A/lOc- 



the agreement becomes worse for smaller values (A/uo c < 
0.1). The discrepancy in A/lu c is typically around 10%, 
but becomes larger for very small A/u> c . 







ou 2 c S(0) 




err. 


(A/oj c ) corr 


err. 




Ising 




Shiba Shiba 




cutoff 


0.05 


324 


1.19 x 10 5 


389 


20% 


0.071 


41% 


0.1 


117 


1.22 x 10 4 


124.6 


6.5% 


0.118 


18% 


0.2 


41.4 


1.35 x 10 3 


41.5 


0.25% 


0.198 


1.0% 


0.5 


8.19 


53.6 


8.26 


0.87% 


0.445 


11.0% 



For the biased system, a gener, 
is given by Lesage and Saleurl 



lization of expression (W 



A| 



lim S(uj) = 

u^o 7r(e 2 + A2) 2 



(42) 



Figure 12 compares this to the numerical results. The 
agreement with Shiba's relation is excellent for small a, 
but only qualitatively right for larger values of a > 1/2. 



IX. CONCLUSIONS 



In the Toulouse case, a =-L/2, the explicit form of the 
spectral function is given aso (note the different normal- 
ization of the spins) 



7rA| 



sH + 4 



, _i / W 

tan - — 

uj \ At 



In 1 



where 



A T 

u 



7T A 

A T = - — . 

4 UJn 




(38) 



(39) 



Figure 13 compares the analytic expression to our numer- 
ical results. The simulation data for 0.1 < A/oj c < 0.2 re- 
produce the analytical result exactly. Outside this range, 
the agreement is still good but no longer exact. 
At zero frequency, expression (|3q) reduces to 



lim S(u>) 

ul—>0 



ttA t ' 



(40) 



Equations (^) and (Eo|) allow us to estimate the cutoff 
corrections to A/uj c (see section ||) through the relation 



A 4 < \ 
— = - (WcXo) 



-1/2 



(41) 



The follo wing table compares the static susceptibility, as 
given by ( |37j ) and calculated from Shiba's relation, and 
lists the corrected values of A/uj c . Shiba's relation is 
satisfied to great accuracy for larger values of A/oj c , but 



We investigated the dynamical behavior of a two-state 
system, coupled to a bosonic environment with linear 
spectral density, in the parameter range 0.05 < A/lu c < 
0.5, a < 0.75 and < e/u; c < 0.2. We found that the 
crossover from oscillatory behavior to incoherent expo- 
nential decay occurs in the vicinity of a = 1/2, which cor- 
responds to the Toulouse limit of the anisotropic Kondo 
model. However, technical difficulties prevent us from ob- 
taining precise answers in the limit of very small A/lu c : 
The strong nearest-neighbor interaction in the Ising pic- 
ture causes the imaginary-time correlation function C(u>) 
to fall off very fast, so that we are effectively left with just 
a small number of data points on which we can base the 
Pade approximation. Simulations on a larger Ising-spin 
system would be necessary to yield stable results. 

The results presented within the quasiparticle approx- 
imation support the conclusion that a c = 1/2 is the cor- 
rect crossover value in the limit A/uj c — * 0, but show that 
a c < 1/2 for finite A/u c . 

For a > 1/3, the spectral function S(ui) does not show 
any peaks on the real axis at finite frequencies, but is cen- 
tered around u> = 0. We have seen that the system nev- 
ertheless exhibits oscillatory behavior beyond that point, 
so that a — 1/3 does not correspond to a crossover value. 

A bias increases the oscillation frequency, in the same 
way as it does for the decoupled system, and favors coher- 
ent oscillations to the extent that a strong enough bias 
can induce a transition from overdamped relaxation to 
weakly damped oscillatory behavior. This is true even 
for rather strong coupling 1/2 < a < 1. 
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